/* 
*** Create figures with establishment level outcomes by share of women at establishment
	// input: estlevel_dataset, estlevel_dataset_full, microregions_municipality_concordance
	// output: sh_mgmt_fem_mgmt_het_by_femxshare, sh_mlext_het_by_femxshare, sh_mlprot_mleave_het_by_femxshare, y_lnwage_oldw_fem_het_by_femxshare, y_ln_tot_het_by_femxshare
*/ 

cap log close
cap log using "$logs/establevel_het_log", replace


*********************************
** CHANGES IN FIRM ENVIRONMENT **
*********************************

* data toggles
local singest = 0 //0: no restriction; 1: must be a single establishment firm
local bothFM = 1 //0: no restriction; 1: must employ both men and women at basline
local covered = 1 //0: no restriction; 1: must be in the geographic coverage of the CBA
local signing = 1 //0: no restriction; 1: must be a signing establishment
local geoglevel = "microregion_mode" //state_mode: for state-yr FEs; microregion_mode: for microregion-yr FEs; stateind: for state-industry-yr FEs; microind: for microregion-industry-yr FEs
local no2011 = 1 //0: include 2011; 1: exclude 2011
local balanced = 0 //0: no restriction; 1: restrict to the balanced panel

	* Load establishment sample (main spells)
	use "$files/estlevel_dataset.dta", clear

	* sample restrictions
	if `singest' {
		keep if bl_singest==1
	}
	if `bothFM' {
		keep if bl_bothFM==1
	}
	if `covered' {
		keep if (bl_covered==1)
	}
	if `signing' {
		keep if (bl_signing==1)
	}
	if `balanced' {
		egen xxx = max(exit), by(fakeid_estab)
		keep if xxx==0
		drop xxx
	}
	
	* time varying fixed effects
	//industry
	gen ind2d_num = floor(ind_mode/1e3)
	//geography
	gen state_mode = floor(municipality_mode/1e4)
	rename municipality_mode municipality
	merge m:1 municipality using "$raw/microregions_municipality_concordance.dta", keepusing(microregion)
	drop if _merge==2
	drop _merge
	//fixed effects
	rename municipality municipality_mode
	rename microregion microregion_mode
	egen indyrFE = group(year ind2d_num)
	egen stateind = group(state_mode ind2d_num)
	egen microind = group(microregion_mode ind2d_num)
	egen geogyrFE = group(year `geoglevel')
	if "`geoglevel'"=="stateind" | "`geoglevel'"=="microind" {
		drop indyrFE
	}

	* for clustering
	egen clustergrp = group(fakeid_estab)

	* event study var 
	gen post=(year>=2015) if year!=.
	gen treatpost=treat*post
	tab year, gen(yearj)
		
	forvalues i = 1/7 {
		gen interj`i'=treat*yearj`i'
	}

	lab var interj1 "2011"
	lab var interj2 "2012"
	lab var interj3 "2013"
	lab var interj4 "2014"
	lab var interj5 "2015"
	lab var interj6 "2016"
	lab var interj7 "2017"

	gen zero=0
	lab var zero "2014"
	replace interj4 = zero

	* Outcome variables
	gen sh_mgmt_fem_mgmt = mgmt_fem/mgmt
	label var sh_mgmt_fem_mgmt 		"Share of female managers (among management)"
	gen sh_mlext_mleav = mlext_mleav/mleav
	label var sh_mlext_mleav			"Share on extended maternity leave (among workers on maternity leave)"
	gen sh_mlprot_mleav = mlprot_mleav/mleav
	label var sh_mlprot_mleav			"Share returning after maternity leave (among workers on maternity leave)"
	
	
	* Heterogeneity
	//share of women at baseline
	gen sharefemale = tot_fem/tot
	gen xxx = sharefemale if year==2014
	egen sharewomen_bl = max(xxx), by(fakeid_estab)
	drop xxx
	egen unique_obs = tag(fakeid_estab)

	/* Baseline % women bins (establishment) */
	forvalues i = 0/4 {
	local j1=`i'*0.2
	local j2=`i'*0.2+0.2
	gen share`i'=(sharewomen_bl>=`j1' & sharewomen_bl<`j2')
	}
	replace share4 = 1 if sharewomen_bl==1
	
	qui count if unique_obs==1
	local denom = r(N)
	forvalues i = 0/4 {
	qui count if unique_obs==1 & share`i'==1
	local numer = r(N)
	di "% in share`i' = " (`numer'/`denom')*100
	}
	
	/* Adjust bins: Baseline % women bins (establishment) */
	forvalues i = 0/4 {
	gen xshare`i' = share`i'
	}
	replace xshare3 = 1 if xshare4==1
	drop xshare4
	
	qui count if unique_obs==1
	local denom = r(N)
	forvalues i = 0/3 {
	qui count if unique_obs==1 & xshare`i'==1
	local numer = r(N)
	di "% in share`i' = " (`numer'/`denom')*100
	}
	
	/*** Heterogeneity plot ***/
	* Plot heterogeneity in treatment effect by share of women at the firm.

	local depvar sh_mgmt_fem_mgmt sh_mlext_mleav 
	local j=1

	foreach var of local depvar {

		if `j'==1 {
			local margin = "share women among managers"
		}
		if `j'==2 {
			local margin = "take extended mat. leave"
		}
		
		forvalues i = 0/3 {
			eststo DDf`i': reghdfe `var' treat post treatpost if xshare`i'==1, absorb(fakeid_estab *FE) cluster(clustergrp)
		}
		
		#delimit ;
		coefplot  (DDf0, aseq("0-19") \ DDf1, aseq("20-39") \ DDf2, aseq("40-59") \ 
					DDf3, aseq("60-100") mcol(sourapple1)ciopts(lpatt(solid)lcol(sourapple1%50)))
					
					, 
		
		vertical keep(treatpost) title("") 
		yline(0, lc(sourapple4) lp(solid)) 
		ytitle("DID estimate: `margin'", size(medlarge)) ylabel(-0.02(0.01)0.03, labsize(medium) gmax gmin) 
		xtitle("Share of women at establishment", size(medlarge)) xlabel(none, labsize(medium))
		legend(off)
		plotregion(fcolor(white)) graphregion(fcolor(white)) graphregion(color(white));
		
		graph export "$figures/`var'_het_by_femxshare.tif", replace;
		
		#delimit cr
		local ++j
	}
	
	local depvar sh_mlprot_mleav
	local j=1

	foreach var of local depvar {

		if `j'==1 {
			local margin = "return from mat. leave"
		}
		
		forvalues i = 0/3 {
			eststo DDf`i': reghdfe `var' treat post treatpost if xshare`i'==1, absorb(fakeid_estab *FE) cluster(clustergrp)
		}
		
		#delimit ;
		coefplot  (DDf0, aseq("0-19") \ DDf1, aseq("20-39") \ DDf2, aseq("40-59") \ 
					DDf3, aseq("60-100") mcol(sourapple1)ciopts(lpatt(solid)lcol(sourapple1%50)))
					
					, 
		
		vertical keep(treatpost) title("") 
		yline(0, lc(sourapple4) lp(solid)) 
		ytitle("DID estimate: `margin'", size(medlarge)) ylabel(-0.02(0.02)0.06, labsize(medium) gmax gmin) 
		xtitle("Share of women at establishment", size(medlarge)) xlabel(none, labsize(medium))
		legend(off)
		plotregion(fcolor(white)) graphregion(fcolor(white)) graphregion(color(white));
		
		graph export "$figures/`var'_het_by_femxshare.tif", replace;
		
		#delimit cr
		local ++j
	}



***********************************************************************
** IMPACT OF CUT REFORM ON ESTABLISHMENT LEVEL OUTCOMES (all spells) **
***********************************************************************

* data toggles
local singest = 0 //0: no restriction; 1: must be a single establishment firm
local bothFM = 1 //0: no restriction; 1: must employ both men and women at basline
local covered = 1 //0: no restriction; 1: must be in the geographic coverage of the CBA
local signing = 1 //0: no restriction; 1: must be a signing establishment
local geoglevel = "microregion_mode" //state_mode: for state-yr FEs; microregion_mode: for microregion-yr FEs; stateind: for state-industry-yr FEs; microind: for microregion-industry-yr FEs
local no2011 = 1 //0: include 2011; 1: exclude 2011
local balanced = 0 //0: no restriction; 1: restrict to the balanced panel

	* Load establishment sample (all spells)
	use "$files/estlevel_dataset_full.dta", clear

	* sample restrictions
	if `singest' {
		keep if bl_singest==1
	}
	if `bothFM' {
		keep if bl_bothFM==1
	}
	if `covered' {
		keep if (bl_covered==1)
	}
	if `signing' {
		keep if (bl_signing==1)
	}
	if `balanced' {
		egen xxx = max(exit), by(fakeid_estab)
		keep if xxx==0
		drop xxx
	}
	
	* time varying fixed effects
	//industry
	gen ind2d_num = floor(ind_mode/1e3)
	//geography
	gen state_mode = floor(municipality_mode/1e4)
	rename municipality_mode municipality
	merge m:1 municipality using "$raw/microregions_municipality_concordance.dta", keepusing(microregion)
	drop if _merge==2
	drop _merge
	//fixed effects
	rename municipality municipality_mode
	rename microregion microregion_mode
	egen indyrFE = group(year ind2d_num)
	egen stateind = group(state_mode ind2d_num)
	egen microind = group(microregion_mode ind2d_num)
	egen geogyrFE = group(year `geoglevel')
	if "`geoglevel'"=="stateind" | "`geoglevel'"=="microind" {
		drop indyrFE
	}

	* for clustering
	egen clustergrp = group(fakeid_estab)

	* event study var 
	gen post=(year>=2015) if year!=.
	gen treatpost=treat*post
	tab year, gen(yearj)
		
	forvalues i = 1/7 {
		gen interj`i'=treat*yearj`i'
	}

	lab var interj1 "2011"
	lab var interj2 "2012"
	lab var interj3 "2013"
	lab var interj4 "2014"
	lab var interj5 "2015"
	lab var interj6 "2016"
	lab var interj7 "2017"

	gen zero=0
	lab var zero "2014"
	replace interj4 = zero

	* Outcomes
	gen y_ln_tot = log(tot)
	label var y_ln_tot 				"Log employment"

	* Heterogeneity
	//share of women at baseline
	gen sharefemale = tot_fem/tot
	gen xxx = sharefemale if year==2014
	egen sharewomen_bl = max(xxx), by(fakeid_estab)
	egen unique_obs = tag(fakeid_estab)
	
	/* Baseline % women bins (establishment) */
	forvalues i = 0/4 {
	local j1=`i'*0.2
	local j2=`i'*0.2+0.2
	gen share`i'=(sharewomen_bl>=`j1' & sharewomen_bl<`j2')
	}
	replace share4 = 1 if sharewomen_bl==1
	
	qui count if unique_obs==1
	local denom = r(N)
	forvalues i = 0/4 {
	qui count if unique_obs==1 & share`i'==1
	local numer = r(N)
	di "% in share`i' = " (`numer'/`denom')*100
	}
	
	/* Adjust bins: Baseline % women bins (establishment) */
	forvalues i = 0/4 {
	gen xshare`i' = share`i'
	}
	replace xshare3 = 1 if xshare4==1
	drop xshare4
	
	qui count if unique_obs==1
	local denom = r(N)
	forvalues i = 0/3 {
	qui count if unique_obs==1 & xshare`i'==1
	local numer = r(N)
	di "% in share`i' = " (`numer'/`denom')*100
	}
	qui unique fakeid_estab
	local denom = `r(unique)'
	forvalues i = 0/3 {
		qui unique fakeid_estab if xshare`i'==1
		local numer = `r(unique)'
		di `numer'/`denom'*100
	}

	
	/*** Heterogeneity plot ***/
	* Plot heterogeneity in treatment effect by share of women at the firm.

	local depvar y_ln_tot 
	local j=1

	foreach var of local depvar {

		if `j'==1 {
			local margin = "log employment"
		}

		forvalues i = 0/3 {
			eststo DDf`i': reghdfe `var' treat post treatpost if xshare`i'==1, absorb(fakeid_estab *FE) cluster(clustergrp)
		}
		
		#delimit ;
		coefplot  (DDf0, aseq("0-19") \ DDf1, aseq("20-39") \ DDf2, aseq("40-59") \ 
					DDf3, aseq("60-100") mcol(sourapple1)ciopts(lpatt(solid)lcol(sourapple1%50)))
					
					, 
		
		vertical keep(treatpost) title("") 
		yline(0, lc(sourapple4) lp(solid)) 
		ytitle("DID estimate: `margin'", size(medlarge)) ylabel(-0.15(0.05)0.15, labsize(medium) gmax gmin) 
		xtitle("Share of women at establishment", size(medlarge)) xlabel(none, labsize(medium))
		legend(off)
		plotregion(fcolor(white)) graphregion(fcolor(white)) graphregion(color(white));
		
		graph export "$figures/`var'_het_by_femxshare.tif", replace;
		
		#delimit cr
		local ++j
	}

	reghdfe y_ln_tot treat post treatpost if xshare0==1|xshare1==1|xshare2==1, absorb(fakeid_estab *FE) cluster(clustergrp)
	egen sumtot = sum(tot), by(year)
	sum sumtot if year==2014
	local denom = r(mean)
	egen sumtotx = sum(tot) if xshare0==1|xshare1==1|xshare2==1, by(treat year)
	sum sumtotx if year==2014 & treat==1
	local numer = r(mean)
	di `numer'/`denom'*100
	


************************************************************************
** IMPACT OF CUT REFORM ON ESTABLISHMENT LEVEL OUTCOMES (main spells) **
************************************************************************

* data toggles
local singest = 0 //0: no restriction; 1: must be a single establishment firm
local bothFM = 1 //0: no restriction; 1: must employ both men and women at basline
local covered = 1 //0: no restriction; 1: must be in the geographic coverage of the CBA
local signing = 1 //0: no restriction; 1: must be a signing establishment
local geoglevel = "microregion_mode" //state_mode: for state-yr FEs; microregion_mode: for microregion-yr FEs; stateind: for state-industry-yr FEs; microind: for microregion-industry-yr FEs
local no2011 = 1 //0: include 2011; 1: exclude 2011
local balanced = 0 //0: no restriction; 1: restrict to the balanced panel


	* Load establishment sample (main spells)
	use "$files/estlevel_dataset.dta", clear

	* sample restrictions
	if `singest' {
		keep if bl_singest==1
	}
	if `bothFM' {
		keep if bl_bothFM==1
	}
	if `covered' {
		keep if (bl_covered==1)
	}
	if `signing' {
		keep if (bl_signing==1)
	}
	if `balanced' {
		egen xxx = max(exit), by(fakeid_estab)
		keep if xxx==0
		drop xxx
	}
	
	* time varying fixed effects
	//industry
	gen ind2d_num = floor(ind_mode/1e3)
	//geography
	gen state_mode = floor(municipality_mode/1e4)
	rename municipality_mode municipality
	merge m:1 municipality using "$raw/microregions_municipality_concordance.dta", keepusing(microregion)
	drop if _merge==2
	drop _merge
	//fixed effects
	rename municipality municipality_mode
	rename microregion microregion_mode
	egen indyrFE = group(year ind2d_num)
	egen stateind = group(state_mode ind2d_num)
	egen microind = group(microregion_mode ind2d_num)
	egen geogyrFE = group(year `geoglevel')
	if "`geoglevel'"=="stateind" | "`geoglevel'"=="microind" {
		drop indyrFE
	}

	* for clustering
	egen clustergrp = group(fakeid_estab)

	* event study var 
	gen post=(year>=2015) if year!=.
	gen treatpost=treat*post
	tab year, gen(yearj)
		
	forvalues i = 1/7 {
		gen interj`i'=treat*yearj`i'
	}

	lab var interj1 "2011"
	lab var interj2 "2012"
	lab var interj3 "2013"
	lab var interj4 "2014"
	lab var interj5 "2015"
	lab var interj6 "2016"
	lab var interj7 "2017"

	gen zero=0
	lab var zero "2014"
	replace interj4 = zero

	*Outcomes
	gen y_lnwage_oldw_fem = lnwage_oldw_fem
	label var y_lnwage_oldw_fem 	"Mean of log female earnings (established workers)"
	gen y_lnwage_oldw_mal = lnwage_oldw_mal
	label var y_lnwage_oldw_mal 	"Mean of log male earnings (established workers)"

	* Heterogeneity
	//share of women at baseline
	gen sharefemale = tot_fem/tot
	gen xxx = sharefemale if year==2014
	egen sharewomen_bl = max(xxx), by(fakeid_estab)
	drop xxx
	egen unique_obs = tag(fakeid_estab)

	/* Baseline % women bins (establishment) */
	forvalues i = 0/4 {
	local j1=`i'*0.2
	local j2=`i'*0.2+0.2
	gen share`i'=(sharewomen_bl>=`j1' & sharewomen_bl<`j2')
	}
	replace share4 = 1 if sharewomen_bl==1
	
	qui count if unique_obs==1
	local denom = r(N)
	forvalues i = 0/4 {
	qui count if unique_obs==1 & share`i'==1
	local numer = r(N)
	di "% in share`i' = " (`numer'/`denom')*100
	}
	
	/* Adjust bins: Baseline % women bins (establishment) */
	forvalues i = 0/4 {
	gen xshare`i' = share`i'
	}
	replace xshare3 = 1 if xshare4==1
	drop xshare4
	
	qui count if unique_obs==1
	local denom = r(N)
	forvalues i = 0/3 {
	qui count if unique_obs==1 & xshare`i'==1
	local numer = r(N)
	di "% in share`i' = " (`numer'/`denom')*100
	}

	/*** Heterogeneity plot ***/
	* Plot heterogeneity in treatment effect by share of women at the firm .

	local depvar y_lnwage_oldw_fem
	local j=1

	foreach var of local depvar {

		if `j'==1 {
			local margin = "mean log wage [women; t>12]"
		}

		forvalues i = 0/3 {
			eststo DDf`i': reghdfe `var' treat post treatpost if xshare`i'==1, absorb(fakeid_estab *FE) cluster(clustergrp)
		}
		
		#delimit ;
		coefplot  (DDf0, aseq("0-19") \ DDf1, aseq("20-39") \ DDf2, aseq("40-59") \ 
					DDf3, aseq("60-100") mcol(sourapple1)ciopts(lpatt(solid)lcol(sourapple1%50)))
					
					, 
		
		vertical keep(treatpost) title("") 
		yline(0, lc(sourapple4) lp(solid)) 
		ytitle("DID estimate: `margin'", size(medlarge)) ylabel(-0.03(0.01)0.03, labsize(medium) gmax gmin) 
		xtitle("Share of women at establishment", size(medlarge)) xlabel(none, labsize(medium))
		legend(off)
		plotregion(fcolor(white)) graphregion(fcolor(white)) graphregion(color(white));
		
		graph export "$figures/`var'_het_by_femxshare.tif", replace;
		
		#delimit cr
		local ++j
	}
	
	reghdfe y_lnwage_oldw_fem treat post treatpost if xshare0==1|xshare1==1|xshare2==1, absorb(fakeid_estab *FE) cluster(clustergrp)
	reghdfe y_lnwage_oldw_mal treat post treatpost if xshare0==1|xshare1==1|xshare2==1, absorb(fakeid_estab *FE) cluster(clustergrp)
	
	
cap log close
